arXiv:1509.07301vl [math.NA] 24 Sep 2015 


THREE-DIMENSIONAL SIMULATION OF 
BIOLOGICAL ION CHANNELS UNDER 
MECHANICAL, THERMAL AND FLUID FORCES 


RICCARDO SACCO^ PAOLO AIROLDI^, AURELIO G. MAURI^ 
AND JOSEPH W. JEROME^ 


Abstract. In this article we address the three-dimensional mod¬ 
eling and simulation of biological ion channels using a continuum- 
based approach. Our multi-physics formulation self-consistently 
combines, to the best of our knowledge for the first time, ion elec¬ 
trodiffusion, channel fluid motion, thermal self-heating and me¬ 
chanical deformation. The resulting system of nonlinearly coupled 
partial differential equations in conservation form is discretized us¬ 
ing the Galerkin Einite Element Method. The validation of the 
proposed computational model is carried out with the simulation 
of a cylindrical voltage operated ion nanochannel with K+ and Na+ 
ions. We first investigate the coupling between electrochemical and 
fluid-dynamical effects. Then, we enrich the modeling picture by 
investigating the influence of a thermal gradient. Einally, we add 
a mechanical stress responsible for channel deformation and inves¬ 
tigate its effect on the functional response of the channel. Results 
show that fluid and thermal fields have no influence in absence 
of mechanical deformation whereas ion distributions and channel 
functional response are significantly modified if mechanical stress 
is included in the model. These predictions agree with biophysi¬ 
cal conjectures on the importance of protein conformation in the 
modulation of channel electrochemical properties. 


1. Introduction and Motivation 

Ion channels connect the intracellular environment and the surround¬ 
ing biological medium, allowing the dynamical exchange of the chem¬ 
ical species that supervise the functions of every cell in the human 
body. Ion channels control electrical signal transmission among ex¬ 
citable cells, notably, muscle cells and/or neurons, and because of their 
fundamental role they have been subject to biophysical and mathemat¬ 
ical investigation for a long time. 

The most commonly used approach to ion channel modeling is based 
on the representation of the cell membrane lipid bilayer as an equiv¬ 
alent electrical circuit including capacitive and conductive elements, 
these latter being in general described by a nonlinear current-voltage 
characteristic. The resulting formulation is constituted by a nonlinear 
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system of ordinary differential eqnations (ODEs) whose best known ex¬ 
ample is the Hodgkin-Hnxley (HH) model [2Sl[2S]- We refer to [HHl [T7] 
for a detailed treatment of ODE-based differential models and their 
application to relevant problems in Computational Biology. 

ODE differential models are widely used because of their limited 
computational cost and ease of implementation. However, their bio¬ 
physical accuracy is limited and often allows to characterize only basic 
properties of the cellular system under investigation, such as the home¬ 
ostatic Nernst potential [22], and/or to reproduce simple experiments 
in Electrophysiology, such as the transmembrane current in voltage 
clamp conditions mm- 

To analyze the fundamental mechanisms that govern ion transport 
across a membrane channel, a higher model complexity is required. 

With this aim, a system of partial differential equations (PDEs) ex¬ 
pressing the balance of mass and of linear momentum for each ion 
species can be adopted to describe the dynamical balance of electro¬ 
diffusion forces acting on chemicals flowing across the cell membrane bi¬ 
layer. The resulting formulation is represented by the Poisson-Nernst- 
Planck (PNP) model for ion electro-diffusion. We refer to [ilIHlHHlIH]. 
for a detailed illustration of the PNP differential system and a discus¬ 
sion of its mathematical and numerical properties. 


Chemical flux 



Eigure 1. Scheme of the general physical problem. 

In the present work we address the three-dimensional (3D) study 
of a generalized version of the PNP model, modified by including the 
effect of electrolyte fluid convective transport, a thermal gradient and 
mechanical forces. The extension of the PNP system accounting for 
fluid velocity is well-known as the velocity-extended PNP model and 
has been the object of extensive mathematical and numerical investi¬ 
gation in [m [32l |33l [52l |12| . In these pages we adopt a wider vision 
of ion transport as schematically depicted in Fig. where D is an open 
bounded domain representing the bio-physical system that includes the 
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ions in motion throngh a physiological fluid under the action of several 
externally applied forces: electrical, thermal, chemical and mechanical. 
Each of these external forces contributes to determining ion flow in the 
medium according to the effect of the following vector fields: 

(1) electric field; 

(2) concentration gradient; 

(3) thermal gradient; 

(4) electrolyte fluid. 

In typical applications of Electrophysiology and Cellular Biology the 
electric field may be the result of the application to the cell of an ex¬ 
ternal voltage drop such as in a patch clamp experiment [131 EZl ED], 
and/or the presence of fixed charged ions embedded in the portion of 
lipid bilayer surrounding a transmembrane channel [231 ESj • The ion 
concentration gradient may be the result of the administration of a sub¬ 
stance such as a drug or a toxin [IHl EH EH EH] , whereas the thermal 
gradient may be the result of the activation of heat-sensitive thermo¬ 
transient receptor potential (thermo-TRP) cation channels exposed to 
capsaicin, a natural ingredient of spicy foods such as red hot chilli pep¬ 
pers and of cold-sensitive TRP channels exposed to menthol [131 EH- 
Ion motion can also be activated by a physiological fluid shear stress 
acting on a mechanosensitive ion channel [30] • Conversely, ion motion 
can determine a modification of the physical properties of the surround¬ 
ing environmnent through the effect of electro-osmotic pressure [311ES] , 
whereas the mechanical action can be determined by the chemicals sig¬ 
nal acting to change the size of the channel or of the vessel in which 
the fluid flows as in the mechanism of neurovascular coupling [211111. 

The 3D simulation of the Generalized PNP equations (for short, 
GenPNP) including electrolyte fluid convective transport, thermal gra¬ 
dient and mechanical forces is a formidable task. The presently avail¬ 
able computer resources and numerical discretization schemes, the Galerkin 
Einite Element Method (GEEM) in our case, open the possibility to 
explore at the cellular scale, to the best of our knowledge for the first 
time, the multi-physics complexity of the phenomena illustrated above. 
The computational structure devised to numerically solve in a fully 
self-consistent manner the GenPNP system is part of the software MP- 
FEMOS (Multi-Physics Finite Element Modeling Oriented Simulator) 
that has been developed by some of the authors [321 HD m. MP- 
FEMOS is a general-purpose modular numerical code based on the use 
of the GEEM implemented in a fully 3D framework through shared li¬ 
braries using an objected-oriented programming language (C++). Sev¬ 
eral kinds of finite elements on simplexes are included for appropriate 
discretization of the various PDFs constituting the problem at hand. 

In particular, the electro-chemical and thermal blocks are solved using 
the exponentially fitted edge averaged (EAFE) piecewise linear finite 
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element scheme studied in [T9] and |58|; the flow of the incompress¬ 
ible electrolyte fluid is solved using the Taylor-Hood finite element pair 
(see |1H|, Chapter 9), whereas the mechanical deformation of the chan¬ 
nel is numerically computed using standard piecewise linear finite ele¬ 
ments. 

An outline of the article is as follows. In sec. |2] we introduce the 
basic notation of the problem at hand, its geometry and partition into 
simplices. In sec. the mathematical model of the GenPNP system 
is discussed in detail. In sec. the solution algorithm employed to 
decouple the various blocks of the full system is described. Sec. is de¬ 
voted to illustrate the weak formulation for each differential subproblem 
whereas sec. [^treats the adopted finite element discretization schemes. 
Sec. [^contains an extensive bouquet of simulations aimed at validating 
the biophysical accuracy of model and computational schemes. Finally, 
sec. [^collects concluding remarks and future research objectives. 

2. Basic Notation, Geometry and Geometrical 
Discretization 

Referring to Fig. [T| in this section we introduce the basic nomen¬ 
clature associated with the multi-physics problem, the computational 
domain in which the problem has to be solved and its geometrical dis¬ 
cretization needed for the subsequent numerical approximation with 
the GFEM. 

2.1. Geometry and notation. Let x and t denote the spatial and 

temporal coordinates, respectively. The simulation domain D C is 
assumed to be an incompressible, viscous and Newtonian fluid mov¬ 
ing with velocity u = u(x, t), in which M > 1 chemical species are 
dissolved. Each chemical has effective charge Zj, i = 1,..., M, and con¬ 
centration (number density) Uj = [m“^]. Cations have Zj > 0, 

anions have Zj < 0 while neutral species have Zi = 0. 

We denote by dQ the boundary of the domain Q and by n the unit 
outward normal vector on dQ. Let U = {n,, ip, u, T, d} denote the set 
of dependent variables of the problem, ip = (p(x, t), T = and 

d = d(x, t) representing the electric potential, temperature and me¬ 
chanical displacement field solutions of the biophysical system under 
investigation. Then, with each variable s £U, we associate a partition 
of the domain boundary dfl into the union of (generally different) sub¬ 
sets. We indicate by Ff) the subset of dfl where a Dirichlet condition 
and by F|,^ the subset where a Neumann condition is applied, in such 
a way that Ff, U F^ = dfl and FI, n F|^ = 0. 

2.2. Geometrical discretization of the computational domain. 

In view of the numerical approximation of a given differential problem, 
the MP-FEMOS computational platform provides a flexible Galerkin 
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Finite Element discretization on a fully unstructured triangulation Th of 
the domain C into tetrahedral elements K, such that [J K = Q. 

KeTh 

For all K G Th, we let hx = diam(E') be the diameter of K, defined 
as the longest edge of K, and indicate hj h = hx the discretiza¬ 
tion parameter. We assume to be a regular partition (cf. |1H] , Chap¬ 
ter 3), i.e., that there exists a positive constant ^ (independent of h) 
such that 

(1) ^ < e VW e % 

PK 

Pk being the diameter of the largest sphere inscribed in K. Condi¬ 
tion Q is equivalent to assuming a bound from below of the minimum 
mesh angle. 


3. Mathematical Model 

In this section we illustrate the multi-physics model for ion elec- 
trodiffnsion in the biological environment schematically represented in 

Fig.[T} 

Ion and fluid motion are mathematically described by a coupled sys¬ 
tems of PDEs including: the velocity-extended Poisson-Nernst-Planck 
(VE-PNP) equation system to account for electro-chemical forces and 
the Stokes equation to account for fluid forces under the assumption 
that nonlinear convective effects can be neglected because of slow fluid 
motion. 

Thermal dispersion into the channel is considered by a heat conduc¬ 
tion equation that allows to determine temperature evolution at all 
X G n. 

To calculate the deformation of the channel wall, the domain Vt is 
regarded as a compressible linear elastic medium undergoing the clas¬ 
sical Navier-Lame theory. At present, no direct coupling is activated 
between stress and deformation fields and the PNP system. This mod¬ 
eling limitation will be removed in a future step of our research. 


3.1. The VE-PNP Poisson-Nernst-Planck system. A detailed 
derivation of the PNP systems can be found, e.g., in [12]. By the 
application of: 1) mass balance for each chemical; 2) linear momen¬ 
tum balance (second law of dynamics); 3) Stokes’ law for viscous drag 
exerted by a fluid on a particle in motion through it; 4) Einstein’s re¬ 
lation; and 5) the law of perfect gases, the following system of PDEs 
is obtained: 
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Vz = 1, 

..., M, solve: 


(2a) 

^ -F div fj = 0 
ot 


(2b) 

Zi VT 

fi — 1 1 Dj^Tli DiTli 

l^zl T 

riiU 



M 


(2c) 

div (e/E) = g ^ ZiUi + qp fixed 

i=l 


(2d) 

E = —Vp. 


In (0), 

the symbol f, denotes the ion particle flux [m ^s ^], pi and 


Di are the ion electrical mobility and diffusivity and T is system tem¬ 
perature. In the Poisson equation ( |^ , E and cp are the electric field 
[Vm“^] and electric potential [V], respectively, while q and e/ denote 
the electron charge and the electrolyte fluid dielectric permittivity, re¬ 
spectively. The quantity qp fixed [Cm“^] mathematically accounts for 
the presence of a fixed charge density due to the lipid membrane bilayer 
and is assumed to be a given function of position only. The diffusion 
coefficient Di and the mobility pi are proportional through the gener¬ 
alized (because T is not constant) Einstein’s relation 

/ON n - 

(3) Di Mi I I 

where ks is Boltzmann’s constant. Initial conditions for ion concentra¬ 
tions are Vi = 1,..., M: 

(4a) nj(x, 0) = n°(x) in 

where the functions are positive given data. The boundary condi¬ 
tions for the VE-PNP system are Vi = 1,..., M: 


(5a) 

p = p 

on 

r£ 

(5b) 

E • n = 0 

on 

^ N 

(5c) 

rii = Ui 

on 

-prii 

^ D 

(5d) 

fi • n = 0 

on 

^ N 


where Tp is the electrostatic potential of the side P^ and n* is a given 
concentration of the chemical species i on side P^h Conditions (5a) 
and (5c) enforce respectively a given voltage and a given concentration 
on the side, while (5b) and (5d) enforce a homogeneous Neumann con¬ 
dition, that corresponds for chemical particles to the conservation of 
mass concentration inside the simulation domain (i.e. particles cannot 
leave the domain). 

The VE-PNP equation system can be regarded as a generalization of 
the classic Drift-Diffusion (DD) model for semiconductor devices [53l 
|39l uni |3T] through the addition of the thermomigration and fluid veloc¬ 
ity terms in the linear momentum balance equation (2b). As a matter 
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of fact, the third term in (2b) is related to the thermal grandient effect 


on ion flow while the last term gives rise to an additional translational 
contribution due to fluid motion and is highlighted by a box. This term 
represents the coupling between electrolyte fluid motion and electrod- 
iffusive ion transport. Henceforth, we use the reference ”PNP system” 
every time the contribution of the velocity and of the thermal gradient 
is neglected in the equation set 0. 

3.2. The Stokes system. The Stokes system to describe the slow 
motion of an incompressible and viscous fluid with a constant density 
Pf [Kg m“^] consists of the following PDEs (see |1H] for a complete 
mathematical and numerical treatment): 

(6a) div u = 0 

, , , (?u , , , 

(6b) Pf—=diy g{n,p) + 

(6c) ^(u,p) = 2/i/e(u) — 

(6d) |(u)=y>) = l(Vu+(Vuf) 

where u is the fluid velocity [ms“^], p is the fluid pressure [Pa], pf 
the fluid shear viscosity [Kg m“^s“^], a the stress tensor [Pa] and e the 
strain rate tensor [s“^]. The symbol 5 is the second-order identity tensor 
of dimension 3 whereas the second-order tensor V^(u) is the symmetric 
gradient of u. Notice that in accordance with the assumption of slow 
fluid motion, the quadratic convective term in the inertial forces has 
been neglected in the momentum balance equation (6b). The boxed 



term at the right-hand side of equation (6b) physically corresponds to 


the electric pressure exerted by the ionic charge on the electrolyte fluid, 
and mathematically represents the coupling between electrodiffusive 
ion transport and electrolyte fluid motion. The initial condition for 
electrolyte fluid velocity is 

(7a) u(x, 0) = u°(x) in 

where the function u° is a given datum, usually set equal to zero. The 
boundary conditions for the Stokes system are: 

(8a) u = g on 

(8b) a;(u,p)n = h on 

where g is the imposed velocity on the side P^ and h = pn is the 


value of the normal stress on P]^. Relation (|8a|) is the wall adher¬ 


ence condition (typically g = 0 ) and (8b) is the application of the 
action-reaction principle for the normal stress an. The presence of the 


two boxed terms in the linearized momentum balance equations (2b) 


and (6b) terms requires the adoption of a suitable iterative procedure 
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to successively solve the whole coupled system (|^- Q as described in 
sec. m 


3.3. Heat equation. Ion and fluid motion, as well as the action of 
mechanical forces, may dissipate energy throughout the channel. The 
evolution of the resulting channel thermal profile T = T(x.,t) in fl is 
obtained by solving the thermal conduction equation (see [HH!) given 
by 

(9a) + div q = Qheat 

(9b) q = —kVT 

where c is the specific heat [w?s~‘^K~^]^ q is the heat flux [ITm“^], k is 
the thermal conductivity [Wm~^K~^] and Qheat is the heat generation 
term [Wm~^]. We notice that convective thermal effects due to fluid 
velocity are neglected in the mathematical definition of the heat flux 
q. Also, we notice that the heat production term Qheat is the result of 
heat source/sink processes such as Joule heating or (possibly nonlinear) 
thermo-chemical reactions due to the administration of drugs as in TRP 
channels. For the purpose of the present work, Qheat is set equal to zero, 
but this limitation will be removed in a subsequent step of the research. 
The initial condition for electrolyte fluid temperature is 

(10a) T(x, 0) = r°(x) in Q 

where the function is a positive given datum. The boundary condi¬ 
tions for the channel heat equation are: 

(lla) T = T on T^ 

(llb) q-n = 0 on T^ 


where T is the imposed temperature on the side T^. Relation (11a) 
represents the effect of an infinite heat source at temperature T = T 
while (11b) is the adiabatic condition. 


3.4. The Navier-Lame approach to mechanical equilibrium. 

The mechanical problem of channel deformation is here described by 
the model of linear isotropic elasticity including thermal expansion and 
represented by the following system expressing the equilibrium condi¬ 
tion and the generalized Hooke’s law: 

(12a) div + fmeo = 0 

(12b) ^(d-. T) = £(|_(d) - |“ ^(T)) + ^ 

(12c) := V,(d) 

where d = [dx,dy, is the displacement [m], is the stress tensor 
[Pa ], fmec is the volumetric force density [N/m ^], is 
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the total strain tensor, is the elastic strain tensor, is the ther- 

—mec —mec 

mal strain tensor, C is the elastic tensor [Pa] and (Tq is the initial stress 
[Pa], In particular, if the thermal properties are homogeneous and 
orthotropic the thermal deformation (strain) is given by the relation 

.. = ai-{T-Tref) 

where ai, a 2 and as [K~^] are the linear thermal expansion coefficients, 
and Tref is the reference temperature [K]. Under standard assumptions 
(symmetric stress and strain, symmetric and positive definite energy 
functionai and coordinate invariance of mechanical response) the rela¬ 
tion between stress and elastic strain becomes 


(13) 


a .. = 2as 

—mec,ij 


el 

mec,ij 


Xs 


el 

mee,kk 


i, j — 1,2,3 


where 


A 


uE 

(1 -h 0(1 - 20 


E 

^“ 2(1 + 0 


are the Lame coefficients, while 0 < < 0.5 is the Poisson ratio and 

U > 0 is the Young module [Pa] describing the material mechanical 
properties. In the mechanical system, at each time level t, the boundary 
conditions are expressed by a given displacement on side or by a 
stress free surface on side P^: 


(14a) d = du on P^ 

(14b) a(d; r)n = 0 on P^ 


where d^ is a given boundary displacement. 


3.5. Fully coupled multi-physical approach. The multi-physics 
description of biological ion nanochannels is hence obtained by com¬ 
bining together all the previous blocks as follows: 


• Mechanical description 

• Ionic transport Q 


(12a) 


• Thermal description (9a) 

• Fluid velocity calculation (j^ 


We define this as the the Fully Coupled Thermal Velocity Extended 
Poisson Nernst Planck Stokes system in presence of mechanical forces 
or, shortly, FC-MF-T-VE-PNP-S model. To the best of our knowledge, 
such a formulation represents the first example of a multi-physics based 
formulation for ion channel modeling and simulation. 


4. Solution Algorithm 

In view of the numerical solution of the EC-ME-T-VE-PNP-S sys¬ 
tem, for a given time Tf > 0 (final simulation time), we introduce the 
time interval It^ := [0,r/] and we divide It^ into a uniform partition 
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of Ntj, > 1 time slabs of width At = Tf/Nr^, such that 

Ufcr^ = Ixf Then, for each k = 0,, NTf, we denote by: 

(15a) Umrc = [dx, dy, dz] 

(15b) Uj.,,p= 

(15c) [u‘,p‘] 

(15d) U*= 


the solution vectors at time level of the Mechanical block and those of 
the Thermal PNP block, of the Stokes block and of the whole coupled 
problem. The solution of the FC-MF-T-VE-PNP-S system using a 
monolithic scheme is a formidable computational task. Thus, the use of 
a decoupled approach is highly desirable. For each discrete time level 
k = 1,..., Ntj, a two-subblock iteration is carried out. In the first step 
of the iteration the mechanical problem is solved for the displacement 
Umec- In the second step of the iteration Umec is used to compute 
the corresponding deformation 6 = 0{^mec) providing a new domain 
configuration in which the remaining thermal, electrochemical and fluid 
blocks of the problem are successively solved. In doing this, temporal 
semidiscretization is performed using different temporal schemes: a) 
one-step methods: Backward Euler (BE) and Trapezoidal Rule (TR); 
b) two-step method: TR-BDE2. We refer to [H], Chapter 11, for a 
detailed discussion and analysis of these methods. 

The flow-chart of the iteration map that is implemented to solve suc¬ 
cessively each block of the FC-ME-T-VE-PNP-S system is illustrated 
in Fig. where k = 0,... NTf — 1 is the temporal loop counter and 
j is the iteration counter in the PNP cycle. In detail, the successive 
solution of the T-VE-PNP-S system proceeds as follows. 


Step 1:: solve successively the PNP system Q and (9a) using 
u = u^, until self-consistency is obtained among n,, (p and T. 
This step is the Thermal-PNP cycle and returns 
Step 2:: solve the Stokes system ([^ with E = — using 

the Uzawa iterative algorithm (see [IHj, Chapter 9) or a direct 
method. This step is the Stokes cycle and returns 


The criterion adopted to monitor the convergence of the solution 
map of Fig. is to stop the algorithm at the first value j* > 0 of the 
iteration counter j such that 


(16) 


7y(i*+l) _ jjU*) 

'^PNP '^PNP 


< toll 


where toll is a prescribed tolerance, and 

p \ 1/2 


(17) 


| W ||2 = 
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New time step 


Figure 2. Flow chart of the solution algorithm. 


is the 2—norm of a vector w G MF. In the numerical experiments we 
set toll = 10“^. The above described solution map is a generalization 
of the classic Gummel decoupled algorithm widely used in contem¬ 
porary semiconductor device simulation tools for the DD model, and 
analyzed in detail in the stationary case in [39l [3T] and in [iO] in the 
time-dependent case. The convergence rate of the generalized Gummel 
algorithm is linear; however, there are two main advantages by using 
such an approach instead of applying, for instance, a fully coupled 
Newton Method. A first advantage is that in the case of the Gummel 
method memory cost is strongly reduced so that also the overall com¬ 
putational effort is much less intense. A second advantage is that the 
Gummel method is more insensitive to the choice of the initial guess 
than Newton’s iteration. This fact is particularly relevant in multi¬ 
dimensional (and, especially, multi-physics) problems, like that under 
consideration, where the construction of an appropriate initial guess is 
in general not an easy task, sometimes even impossible. 
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5. Weak formulation and well-posedeness analysis 


In this section we discuss the weak formulation and the well-posedeness 
of each differential subblock in the solution algorithm of sec. i In 
sec. 


5.1 we study the T-VE-PNP system, the Stokes system in sec. |5.2 


and the Navier Lame formulation for the mechanical displacement in 
sec. 5.3 The analysis of the Poisson equation (2c) and heat equa¬ 
tion is standard and can be found, e.g., in [271 SH]. For ease 
of notation, we set := ^(x, t^), k = 0, ...,W/; for any function 
^ = ^(x, t). To avoid technicalities in the subsequent analysis we as¬ 
sume, without loss of generality, that the Dirichlet data n* (T-VE-PNP 
system), g (Stokes system) and d/j (mechanical system) are equal to 
zero. Moreover, we restrict the study to the case where the BE method 
is used for time discretization as the analysis proceeds in a similar man¬ 
ner if the TR or the TR-BDF2 scheme is employed. 


5.1. T-VE-PNP system. After time discretization, by using the Ein¬ 
stein relation ([^ in (2b) and the boundary conditions (2a) and (2b), 
the T-VE-PNP system takes the following form: Vfc = 0,..., V — 1 and 
for all j > 0 until convergence of the PNP cycle of Fig. solve: 


+ div fo+i = 0 


(18a) 

7+1 r 

n</ — 

At 
(18b) 

= —Di 


kuT^ 




(18c) 

ry' ■ n = 0 
(18d) 

= 0 


in 


in 


on r)(; 


on 


where is a given function computed by solving the Poisson equa¬ 
tion (2c), is the solution of the heat equation 0 and Dl is the 
diffusion coefficient of species i computed using (j^ with T = TK 


The equation pair (18a)-(18b) is of the admissible form (8.1) of |2()j . 
However, it is immediate to check that the second a priori estimate in 
formula (8.6) of [20] is, in general, not satisfied preventing a theoretical 
study of the well posedness of the equation system (18). To overcome 


this impasse, and also in view of the finite element discretization of (18), 


we replace the expression (18b) with a suitable approximation. To this 
purpose, for any function g : 7/i —t M such that q G W^{K) for all 
K E Thi denote the restriction of q to an element K E Th hj q\K 
and define the piecewise gradient operator : 7/i —)■ (M)^ as the 
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vector-valued function such that 

{^ThO) \k := V (glx) 


'iK e Th. 


Finally, we associate with the function g|x its /larmonic average H^iq) 
defined as 


:= 


Ik 1-' dK 


IA1 


-1 


VK e %. 


Then, the approximate form of the flux (18b) that is used in this article 
is 

(19) 


ff+i = -D^ ( 


qzi 


J+1 

-ni 


Vif' 




j+i 


VTA 


- 




The above definition corresponds to having replaced in the electrical 
drift term of the flux the local environmental temperature with its har¬ 
monic average. It is well known (see H) that the use of the harmonic 
average for rapidly varying functions provides much more accurate re¬ 
sults than the usnal integral average. It is expected that as the mesh 
size /i ^ 0+ also the modeling error introduced by using (19) instead 
of (18b) decreases accordingly. 

We now show that the equation system (18) in which the flux (18b) 
is replaced by ( [I^ is uniquely solvable. Let us introduce the following 
Hilbert space (see [3H! Chapter 1) 

(20) W* = //p",(5d) = {ve H\n) : = 0}. 


Multiplying (18a) by a test function v G integrating over fl and 
enforcing the boundary conditions (18d) and (18c), we obtain the weak 
formnlation of the T-VE-PNP system: V/c = 0,.. .,N — 1 and for all 
j > 0 nntil convergence of the PNP cycle of Fig. find G 
such that 


(21) ^ , v) + o"* (n^+^, v) = (u) 


Vu G C”* 


where: 


(n^+\u) = 


n: 


j+i 


a 




V dfl 


,v)= I • Vu dn + 


qDjzi 


J + 1^ 

-ni 


■ Vv dn 


+ 


■ ■ \7T^ 

• Vu dO, 


Jo, 

r“iv) = I 


Ti 
n^v dVt. 


kBHr,(Ti) 

f -Vv dfl 


We notice that letting At —)■ -|-oo allows ns to recover the steady-state 
problem. Using a standard technique in the analytical and compnta- 
tional study of the DD model for semiconductors, we introduce the 
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change of dependent variable 

( 22 ) hi = nie~^ 


where : 71i —?• M is a dimensionless potential such that G 

for each K eTh, defined as 


qz 


+ ln 


T\k 


'^1'^ - - LSJmO' 

Tref being a reference temperature (typically T^ef = 
ing ( |^ into the expression for the flux ( |18b| ) yields 

(24) 


^ ref / . 




'iK^Th 

300K). Substitut- 


The equivalent form of the flux shows that is composed by the 
sum of two terms, one representing the Pick law of diffusion of the 
novel species with the modified diffusion coefficient Di^e^\ the 
other representing passive convection of driven by the electrolyte 
fluid velocity. 


Remark 5.1. We notice that in the isothermal case the potential 4/ 
coincides with the electric potential ip and the change of variable { 22) 
is nothing but the standard use of the so-called Slotboom variable very 
common in the study and numerical solution of the DD model (see [SSI 

ESI ED;. 


We also make the hypothesis that, Vj > 0 until convergence. 


0 < < Di{x,t) < Df 


almost everywhere in Q. 


Now, replacing (p^ and (p3|) in (21), we obtain the weak formulation 


of the T-VE-PNP system: \/k = 0,..., — 1 and for all j > 0 until 

convergence of the PNP cycle of Fig. 
find hj~^^ G F”'® such that 


(25) 

with 


(nf, v) = F”® (p) Vp G F”® 


1 


6”®(A/',p) = ^ / We^' p d^+ f D^e'^WA/' • Vp dn 

(26) 

- / • Vp dVL VW, p G F’^b 

Jn 


Theorem 5.1. The weak problem (25) is well posed and its solution 
satisfies the following stability estimate 


(27) 


\ni 


1 


< e“ 


HWq 




mm 


Proof. To prove the well-posedness of we apply the Lax-Milgram 
Lemma (see [Hj, Chapter 5). We have: 
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• is a continuous linear functional 

(28a) |f”<(«)| = i 

• •) is a continuous bilinear form 


/ 


n^v dfl 


/ 1 ,, fc, 

-At" * 


lolbllo yveV^-, 


(28b) 

+ Me'^’”“3|AA|M|Vt;||o 

^ ^ ^ | | | | 1 | 1| | 1 

yAf,ve 


• •) is a coercive bilinear form. Given a function / G L^(G) 

we denote by f~^ its positive part. Moreover, we indicate by H 
the Heaviside function such that H{0) = 0. We have 


(28c) 


b-i(N,N) > —e' 


2 + D"e'^'"''*||VW||2 


— e 




IQ 


{Nu’^-VNy dn 


where || • ||o is the L^-norm, || • ||i = ( II ' Ho + I A(‘)llo )Ds the 
i7^-norm. Set / = (WW • VW)^ and notice that / G L^(f2) 
by the Sobolev and Holder inequalities. Using the fact that 
div = 0, a direct calculation yields 


div 




f , 


from which, using the fact that u'' = 0 on we obtain 

[ {Nu'^ • VA^)^ dn = [ f dn = [ div (-N^u’^Hif)] dQ = 0. 
Ja da ka \ 2 / 


Therefore, (28c) yields 


(28d) b^^{N,N)>e 




1 


nim<! Df 


\/N G Wb 


Having verified that properties (28) are satisfied, the use of the Lax- 
Milgram Lemma allows us to conclude that the weak problem (25) 
is uniquely solvable. The stability estimate (27) immediately follows 
from ([28a) and (|28d|). □ 
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5.2. Stokes system. Applying the time discretization to (|^ and en¬ 
forcing the boundary conditions (|^ , we arrive to the following station¬ 
ary generalized Stokes problem: 


(29) 


Pf 


div (2/j/|(u^+i) 


(u/o+i _ 

At 

div = 0 

^(u^+i) = i(Vu^+^ + (Vu^+^)^) 

fc+i ^ 0 


U 




in Q 
in fl 


on 

on T'^ 


where := and being given func¬ 

tions computed by solving the T-VE-PNP system (18) and the Poisson 
equation (2c), respectively. 

We set := [i^Q pu and := L‘^{Q). Then, we multiply 

the first and the second equations of (29) by a test function v G 
and q & QP, respectively, and integrate over fl, to obtain the following 
abstract problem: 

V > 0, A: = 0,..., A — 1, find G and G such that: 


(30a) 

(30b) 

where: 


a“(u^+\ v) + = F'^(v) 

6"’P(u^+\g) = 0 


Vv G 
MqeQP 


a“(u^+\v) = 


U^+l.V 


At 


dO + / 2 p^(u'^+^) : ^(v) dn 
Jn 


^u,p(v^pfc+l) = _ f /+ldiv 

Jq 

F“(v) = 


IQ 

f fc+l . 


dQ + 


V dfl 


h 

• V 


pf 


At 


dO. + 


h ■ V 

Pf 


dV. 


Proposition 5.1. The weak problem (30) is well posed and its solution 
satisfies the following stability estimate 


(31) 


U' 


A:+l| 




. <C„, ||/+‘|la, <C, 


where ||v||yu = ||e(v)||o and HqIIqp = HqUo for all v G P“ and q G Q^, 
while Cu and Cp are positive eonstants depending on problem data, on 
the eoercivity eonstant of the bilinear form a“(-, •) on ker P“ and on the 
inf-sup constant of the bilinear form 6“’^(-, •). 


Proof. Well-posedness of (30) can be proved using the abstract theory 
for saddle-point problems developed in [IH], Section 7.4. The stability 
estimate (31) is a consequence of the application of Theorem 7.4.1 
of [IH]. □ 
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5.3. The mechanical problem. To proceed with the weak formula¬ 
tion of the Navier-Lame model in presence of a thermal field we set 


(32) 




multiply (12a) by a test function v G and integrate over fl. Ap¬ 


plying the boundary conditions (14), the weak formulation of the me¬ 
chanical problem is: 
find G such that: 


(33) 

where: 


^^k+i ^ v) = (v) Vv G 


a'*(d'‘+',v) = ^ ft,.„(d‘+') : |(v) dn 


F^(v) = f (T^) : e(v) dQ + [ div do • v dfl -|- f f^ec • v dQ. 

Jn “ Jn ~ Ja 


Proposition 5.2. The weak problem (33) is well posed and its solution 
satisfies the following stability estimate 


(34) 


IM 


fc+l I 


yd < Cd 


where ||v||yd = ||e(v)||o and Cd is a positive constant depending on 
problem data and on the coercivity constant of the bilinear form •). 


Proof. The proof follows from the application of Lax-Milgram Lemma 
and Korn’s inequality (see, e.g., |H], Thm. 3.4). □ 


6. Galerkin Finite Element Approximation 


In this section we discuss the Galerkin Finite Element approximation 
of each weak problem introduced and analyzed in sec. In sec. |6.1| we 
study the T-VE-PNP system, the Stokes system in sec. |6.2[ and the 


Navier Lame formulation for the mechanical displacement in sec. 6.3 


Throughout this section, for r > 1, we denote by Pr(K) the space of 
algebraic polynomials of degree < r defined on the element K £ Th- 


6.1. T-VE-PNP system. Let us define the following finite-dimensional 
subspace of ( |^ 

w* = {Vh G G°(n) : Vh\K e Pi(K) VK G %, = 0} 


such that Nh = dim (WO- The Galerkin Finite Element approximation 
of (25) in matrix form is the linear algebraic system 
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where is the generalized stiffness matrix, M*** is the 

mass matrix, is the stiffness matrix and is the load vector, 
given by: 


(36a) 

(36b) 

(36c) 


Ar ^ ^ 

F"* = F”'($i)- 


In (|36|), is the Lagrangian set of basis fnnctions for whereas 


6)(*(-, •) is a suitable modification of the bilinear form (26). The solution 


vector of the linear system (35) contains the degrees of freedom 
(dofs) of the approximation of the unknown ion concentration n* 
at t im e 

Three computational issues are worth noting in the stable and accu¬ 


rate solution of system (35). The first issue concerns the mass matrix 
M"*. We compute its entries using the 3D trapezoidal rule which corre¬ 
sponds to ” lumping” diagonalization of M"*. The second issue concerns 
the numerical quadrature rule used to compute the entries 6/j(d>j,$j) 
of the stiffness matrix As a matter of fact, if the convective term 
of the flow fn; dominates the diffusive contribution, the standard fi¬ 
nite element method with piecewise linear polynomials may lead to 
instability, with spurious oscillations in the numerical results. An es¬ 
tablished remedy to this drawback is the EAFE (Edge Averaged Einite 
Element) approach studied in 2D in [JH] and in 3D in [SH!- This latter 
is a multidimensional generalization of the celebrated ID Scharfetter- 
Gummel method [51], and is obtained by replacing along each edge of 


the element K the diffusion coefficient in (26) with its har¬ 


monic average. The third issue concerns the choice of the dependent 


variable in the computer implementation of the linear system (35). As 


a matter of fact, the change of variable (22) gives rise to a symmet¬ 
ric positive definite stiffness matrix B'^\ However, the entries of this 
matrix may be impossible to compute because of overflow exceptions 
due to the dynamic range required by the evaluation of the quantity 
Thus, it is necessary to go back to the original dependent variable 


o’!- 


n by applying the inverse of the transformation (22) at each node of 
Th- This is equivalent to multiplying each column I of H*** by for 
I = 1,..., Ai/j (cf. [Ill [H]). The novel linear algebraic system for the 
variable n can be written in matrix form as 


(37) 

where 




k^l _ 


= X^Miag(e 




Solving (37) yields an accurate and robust treatment of the mass bal¬ 


ance equations for rii especially in the presence of high convective terms 
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(E and/or u), preventing the occnrrence of spurious unphysical oscil¬ 
lations in the computed ion concentrations under suitable conditions 
on the regularity of Th- 

Proposition 6.1 (Discrete maximum principle). Assume u = 0 and 
that the triangulation Th satisfies Lemma 2.1 of [HH] • Then, is a di¬ 
agonally dominant by columns M-matrix. This implies that system (j^ 
is uniquely solvable [HS] and that if¥^^^ > 0 then > 0 (Discrete 
Maximum Principle, DMP). 

The previous result shows the beneficial properties of the extension 
of the EAFE method to the numerical approximation of the PNP-T 
model. The analogue of Prop. |6T] in the case of the VE-PNP-T model 
(where the electrolyte fluid velocity is non vanishing) is at moment 
an open issue that will be the object of a next step of the current 
research. We limit ourselves to the observation that in all the numerical 
experiments reported in sec. no spurious oscillations ever appeared 
in the computed ion concentrations. This is to be ascribed to the fact 
that u is a small perturbation to thermo-electrochemical flow so that, 
in practice, the method continues to satisfy a DMP. 


6.2. Stokes system. In sec. |6.2.1| we describe the two finite element 
schemes that are used in MP-FEMOS for the numerical discretization 
of the saddle-point problem (30) whereas in sec. 6.2.2 we illustrate and 
analyze the algorithms that are implemented for the efficient solution 
of the linear algebraic system. 


6.2.1. Finite element schemes. The Galerkin finite element approxima¬ 
tion of the generalized Stokes problem (jS^ is a delicate issue because 
if the finite element spaces for the electrolyte fluid velocity and pressure 
do not satisfy a compatibility condition (the so-called inf-sup or LBB 
condition [10]) then uniqueness of the discrete solution fails to hold 
and spurious pressure modes may arise affecting the stability of the 
formulation. To overcome this difficulty two general approaches can be 
adopted. The first approach consists in selecting a finite element pair 
satisfying the inf-sup condition 1) in a manner that is independent of 
the discretization parameter h; and 2) with an optimal convergence or¬ 
der in the graph norm with respect to the V xQ topology. The second 
approach consists in modifying the saddle point problem (30) by the 
introduction of a stabilization term whose weight is strong enough to 
circumvent the inf-sup condition. 

The first approach is pursued in MP-FEMOS by implementing the 
so-called Taylor-Hood pair (see [Hj Chapter 9 for a detailed description 
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and analysis). This amounts to setting: 


(38a) 


= {v;, G [C%n)f : G [P2(i^)]' yK G r,,v,|ru = 0} 

(38b) 

Qr = Uh e c%n) : QhlK e Pi(i^) yK G Th} 

and to solving the following sequence of saddle-point problems: 

V > 0, /c = 0,— 1, find G 14^^ and G such that: 

,k+l 


(38e) a"(u*+',Vj) + (,"-'’K.p„) = F"(vi) Vv^ e V™ 

(38d) (>“-f(u‘+>,®)=0 V*eQ™. 

The second approch is pursued in MP-FEMOS by implementing the 
so-called Hughes-Franca-Balestra stabilization (see [2H])- This amounts 
to setting: 

(39a) 

i/HFB 

(39b) 


= K e [C'°(b!)]3 : v,|k e [Pi(iF)]3 ViF G = 0} 


gr " = {Qh G C%Q) : qh\K G Pi(iF) ViF G %} 


and to solving the following sequence of stabilized saddle-point prob 
lems: 

y = -I, find G and G Q 

that: 


such 


(39c) 

(39d) 

where: 


^hfb 

h 


a"(uy‘.vO + 6"'!'(v,„Pj) = F"(vJ 


(®) := 


hi 


K 


fK 


- pAu^+ 1 
At /i 


Vv, G Vf 
Vg. G 

Vph - f 1 • Vqh 


f =-h 

Pf 

The quantity (5 is a positive parameter to be properly selected upon 
noting that a too large value of 5 does not eliminate the spurious modes 
of the pressure while a too small value of d yields a poor approximation 
for the pressure field near the domain boundary (cf. [IH], Chapter 9). 


6.2.2. Analysis and linear solvers for the generalized Stokes system. In 
this section we assume that 1) the time advancing scheme has second- 
order accuracy (as in the case of the Crank-Nicolson scheme or the 
TR-BDF2 scheme); and 2) the exact solution is appropriately smooth. 
The quantity C denotes a positive constant independent of h whose 
value is not the same at each occurrence. 
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The Taylor-Hood FE solution U' 


TH 


following optimal error estimate (see [IHj, Chapter 9) 
(40) ||u^+^-u' 


X G satisfies the 


TH\ 




TH\ 


In spite of the fact that ( |40| ) shows that the TH method is second- 
order accurate, the linear system associated with the saddle-point FE 


problem (38c)- (38d) has a very large size in 3D computations, so that 


it is a good practice to use an iterative method for its solution. In MP- 
FEMOS, the adopted approach is the Uzawa method (see |18] Chapter 
9) with preconditioner P given by the mass matrix scaled by the fluid 
viscosity u and acceleration parameter p such that p > l/{2u). In 
the case where convergence of Uzawa’s iteration is still to slow, an 
alternative is given by the MUltifrontal Massively Parallel sparse direct 
Solver (MUMPS) developed at INRIA 


The HER FE solution uf G 


■f/HFB 


X 


G satisfies the 


following error estimate (see [Hj, Chapter 9) 


(41) 


u' 


Ai+1 


— U 


HFB\ 


\V 


+1 - 


Ph I 


Iq < Ch. 


The above result is suboptimal for the pressure variable. However, the 
HFB stabilization, compared to the TH approach, has the advantage 
of allowing the use of equal-order interpolation for both the unknowns 
of the Stokes system. This reflects into a simpler coding, a much re¬ 
duced memory size and, more importantly, into a substantially lower 
computational effort. The disadvantage lies principally into the need 
of a proper selection of the stabilization parameter 5. 

A thorough comparison between the performance of the two solution 
methods considered in the MP-FEMOS software is addressed in sec. [71 


6.3. The mechanical problem. Let us define the following finite¬ 
dimensional subspace of (32) 

= K e [^“(n)]^ : e [Pi(/^)]' ViF G r,, v,|rd = 0 }. 


The Galerkin approximation of 
Find G such that: 

( 42 ) 


is: 


a''(dy',vj)=F»(vk) Vv»el4» 


Problem (42) is equivalent to the solution of the linear algebraic system 

(43) = L‘^ 

where: 


4 = 

L? = F«(4.j) 


b i — b •••) Nt 
i = 1 ,..., Nt, 


Nt being the number of dofs for each component of the displacement 
vector d. The stiffness matrix is symmetric and positive definite 
from which it follows that (43) is uniquely solvable. No ill-conditioning 
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due to material compressibility affects the linear elastic mechanical 
problem because the material Poisson ratio v is far from the limit value 
0.5 (cf. Tab. [^. Using standard arguments (see [2Zli) the following 
error estimate can be proved 

||u - u/i||yd < C^h 

where C^. is a positive constant independent of h and depending only 
on the mesh aspect ratio 7 , on problem data and on the exact solution 


of (33). The numerical solution of the linear algebraic system (43) is 


efficiently carried out in the MP-FEMOS platform through the use of 
the Bi-CG iterative method (see 133 . Chapter 4 for a discussion of the 
algorithm and its convergence analysis). 

7. Simulation Results 

In this section we present the results of the reciprocal interaction 
among the thermal and mechanical aspects and the electrodiffusion of 
ions in a nanoscale biological channel. The considered case study is the 
same channel already investigated in two spatial dimensions in [ 21 ] and, 
more recently, in three spatial dimensions in [I] . The biological setting 
is constituted by a voltage operated channel in which K"*" (potassium) 
and Na'*' (sodium) ions are taken into account. In the hrst step of our 
analysis a numerical simulation of the VE-PNP model is performed 
to investigate the coupling effects between the PNP and the Stokes 
systems in a 3D cylindrical geometry used as a natural extension of a 
2 D domain as already demonstrated in [T]. In the second step of our 
analysis we carry out the study of the influence of the thermal gradient 
in the VE-PNP model. In the third step of our analysis a mechanical 
stress (with different strengths) is applied at the center of the cylinder 
causing a channel deformation. Under such a condition we investigate 
ion electrodiffusion using the VE-PNP model to highlight the variation 
in ion flow induced by channel restriction. Eor completeness of infor¬ 
mation, Tab. [^summarizes the values of all model parameters and also 
the initial and boundary conditions for the VE-PNP system. 

In Figs. 3(a) and |3(b)| the mesh structures used in the numerical 
simulations are reported. The typical number of simplices is in the 
range of 50000 elements. Channel length is set at lOnm while channel 
diameter is equal to 2 nm. Channel terminals are restricted to the 
internal cylinder of the lateral vertical surfaces: the outside structure 
represents the biological region that is assigned to the constriction or 
to the dilation of the channel itself. The Z-axis is aligned with the 
symmetry axis of the structures and is oriented from Side A towards 
SideB terminals. Cell interior is located at SideA while cell exterior 
ambient is located at SideB. 


7.1. Numerical results for the VE-PNP model. In this case study 
we investigate the impact of the velocity term in the VE-PNP model at 
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Parameter 

value and units 

Zk+ 

+1 


+ 1 

K~^\sideA 

2.41 • 10^° cm“^ 

NcL'^ Is'zdeA 

3.01 • 10^® cm“^ 

K~^\sideB 

1.2 • 10 ^® cm“^ 

Na+\sideB 

2.65 • 10^0 cm-3 

I^K+ 

7.2- 10-^cm^V-^s-^ 


5.2 • 10 “^ cm^V“^s“^ 


1.2 • 10 ^^ cm“^ 

^init 

^Na+ 

3.01 • 10^*^ cm“^ 

^ SideA 

0.02 V in sec. 17.11 and 17.21 0.2 V in sec. 17.31 

^IsideB 

OV 

P\ SideA 

0 Ncm“^ 

P\sideB 

100 Ncm“^ 

flf 

10“^ Nscm“^ 

Pf 

10“^ Kgcm“^ 

E in Qi and 

1.5 • 10^ Ncm“^ 

E in Q 2 

1.5 • 10^ Ncm“^ 

E in Cylinder 

7 • 10^ Ncm“^ 

u in r2i and ^3 

0.3 

V in Vt 2 

0.4 

V in Cylinder 

0.2 

T 

293.75 K 


Table 1. VE-PNP model parameters in presence of 
thermal and mechanical forces and boundary conditions 
for Poisson and electrodiffusion equations. 



Figure 3. Mesh used in the simulations: (a) channel 
with no deformation; (b) channel with deformation due 
to the application of an internal stress in Q 2 - 




































MCCARDO SACCO\ PAOLO AIROLDl\ AURELIO G. MAURl\ AND JOSEPH W. JEROME^ 

a constant room temperature. Referring to ([^ and (|^, and according 
to Tab. Dirichlet boundary conditions for ionic concentrations and 
electrostatic potential are enforced on SideA (cell-inside) and SideB 
(cell-outside) surfaces while homogeneous Neumann conditions hold 
elsewhere. In the Stokes problem i§ Neumann boundary conditions 
are applied on SideA and SideB, with a pressure drop of 800Ncm“^ 
whereas homogeneous Dirichlet conditions for fluid velocity hold else¬ 
where (adherence condition). The initial values are set at and 

(7^*+ for the ionic concentrations and equal to zero for the electrolyte 
fluid velocity. 




(a) Electrostatic potential. 


(b) Electric field. 


Figure 4. Comparison between VE-PNP and PNP 
models of the electrostatic potential (a) and electric field 
(b). ID cut is along the channel axis at t = 50ns. 


Fig. 4(a) and fig. |4(b)| show at time t = 50ns the profile of electro¬ 
static potential and of the 2 : component of the electric field obtained 
along the cylinder symmetry axis for both VE-PNP and PNP-only 
models. The onset of pronounced bumps in the electrostatic potential 
is responsible for a nonuniform contribution of the electric field to ion 
electrodiffusion. As a matter of fact, close to the interior part of the 
channel drift forces are directed towards cell interior while at the end 
of the channel the drift forces point outward of the cell. The difference 
between the prediction of the VE-PNP and PNP models is not so rele¬ 
vant in this case, as confirmed by the fact that the small variations in 
the potential plot completely disappear in the electric field plot. 

Pig.j^shows the profile of K"*" and Na"*" ion concentrations. Since the 
^-component of the electrolyte fluid velocity is positive the K"*" species 
is more diffused towards the positive 2 ;-axis in the case of the VE-PNP 
model. Conversely, in the case of Na"*" the electrolyte fluid flow causes 
a reduction of diffusion towards cell interior. The asymmetry in the 
ionic concentration profiles as shown in fig. is to be ascribed to the 
cumulative effect of drift and diffusion forces: in fact, close to SideA 
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Coordinate z [nm] 


Figure 5. Comparison between PNP and VE-PNP 
models of the K+ and Na"*" ionic concentrations. ID cut 
is along the channel axis at t = 50ns. 

and SideB drift is opposing to diffusion but K"*" ions have a higher 
diffusion gradient and lower drift force. In fig. a vector field plot of 
the electrolyte fluid velocity is shown at time t = 50ns: as anticipated, 
velocity is directed towards the positive .s-axis due to the imposed 
boundary conditions. 





Velocity: Z-Component (um s-1) 

0 1e+5 2e+5 


209641.4 


-0.01025 


Figure 6 . VE-PNP model: 3D view of electrolyte fluid 
velocity vector field inside the channel at t = 50ns. 


7.1.1. Validation of the Hughes-Franca-Balestra formulation. As already 
pointed out in se c. |6 .2| the use of a stabilized finite element space in the 
discretization of Q considerably reduces the number of dofs and, as a 
consequence, the computational cost. The HEB methodology, however, 
requires an appropriate choice of the stabilization parameter S. This 
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issue is documented in the example below. Fig. 7(a) and fig. 7(b)| show 
the velocity and the pressure profiles in a case study on a cylinder of 
lOOnm length and 25nm of diameter where on SideA the 2 :-component 
of the velocity is set equal to 0.1/ims“^, on SideB the Neumann condi¬ 
tion aja. = 2 • 10“^nPa is applied while a no-stress condition is applied 
on the lateral surface. The numerical solutions obtained using MP- 
FEMOS with the TH or HFB pairs are indistinguishable from those 
computed by the commercial software COMSOL 




(a) (b) 

Figure 7. Comparison between the different discretiza¬ 
tion spaces for eq. Q and the result obtained using the 
commercial software COMSOL [Hj: (a) velocity along 
(b) fluid pressure. 


7.2. Effects of a thermal gradient in the VE-PNP model. In the 

context of the numerical simulation of ion migration in phase change 
materials for electronics applications sa we have introduced the effect 
of a thermal gradient in the physical system as a further driving force 
in the basic electrodiffusion model of ions (2a). In this section we aim 
to understand the role of such phenomenon in the context of a bio¬ 
logical environment. To verify this issue we perform a computational 
test in which we have deliberately chosen a non-physiological value of 
the temperature located outside the cell {SideB terminal). The ther¬ 
mal gradient is consequent result of the following boundary conditions 
in (lllal): 


T = 293.75 K 

on Bot 

T = 343.75 K 

on Top 

T = 373.75 K 

on Top 

VT • n = 0 

on r)0. 


We named the previous conditions as easel and case2 respectively 
where we have neglected any heat generation phenomena in (9a). The 
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thermal profiles found inside the ionic channel after the simulation 
are complicated by the fact that the Dirichlet boundary conditions are 
applied only on the channel terminals and not on the whole T as shown 
in Figs. |8(a)| and |8(b) , 



(a) 



(b) 


Figure 8. Channel thermal profiles resulting from (9a): 


(a) 3D view of easel; (b) ID-cut along the geometrical 
symmetry axis for both easel and case2. 


Figs. 9(a) and |9(b)| show the effect of the thermal gradient on K+ 
and Na’*' concentration profiles: ID-cut are at the final time Tf = 
50ns along the geometrical symmetry axis. Black color refers to the 
constant temperature case T = 300K, the red color to easel and the 
blue to case2. As the thermal gradient increases, the ionic species are 
drifted towards the SideA surface, as prescribed by the negative sign 
in eq. (9a): Na"*" is more diffused while K"*" is less. 

The electrostatic potential, shown in fig. |10(a)| is increased in the 
left part of the channel due to the higher ion concentrations in that 
region when the thermal gradient is not null. For the same reason the 
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(a) 



(b) 


Figure 9. K+ (a) and Na+ (b) concentration profiles: 
red color constant temperature at T = 300K, blue color 
easel and black color case2. 


infiuence of the ionic pressure onto the fluid is more pronounced as 


shown in fig. 10(b) 


From the computational analysis of this section we can conclude that 
the contribution of a thermal driving force to ionic flow is not relevant in 
biological systems unlike the case of semiconductor materials illustrated 
in 


7.3. Numerical results for the VE-PNP model considering me¬ 
chanical stress. In this section we demonstrate the effect on ion elec¬ 
trodiffusion due to the action of external mechanical forces able to 
modify the size of the biological channel. This aspect of cellular bi¬ 
ology is fundamental in many control mechanisms of a living system, 
for example, heart beat, muscle activity or cerebrovascular autoregula¬ 
tion (see [35] for details and further reference). As discussed in sec. 
the solution of the mechanical problem (12a) is used to determine the 


finite element mesh transformation that yields the new domain geom¬ 
etry in which ([^ and (|^ are solved with a split domain algorithm. 
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(a) 



(b) 


Figure 10. Electrostatic potential (a) and . 2 - 
component of the fluid velocity (b) profiles: red 
color constant temperature at T = 300K, blue color 
easel and black color case2. 


The following boundary conditions are applied to solve (12a) where 
Tpf is the lateral surface of the whole parallelepiped, fli, VL 2 and fla are 
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three separate regions with different mechanical properties as reported 
in Tab. [7l 


(44a) 

dn = [0,0,0]^/im 

(44b) 

fmec = [0, 0,0]^ Ncm“^ 

(44c) 

^ = [ 7 , 7 , 0 , 0 , 0 , 0 ]'^ Ncm"^ 
(44d) 

7 = -1 • 10 ^ 

(44e) 

7 = -2 • 10 ^. 


on SideA U SideB UFlUTbUTt 
in fli U U Q 3 
in fl 2 


We denote condition (44d) as deformation 1 (de/.l) while (44e) as 
deformation 2 {def. 2 ): the case with no deformation is indicated with 
no — def.. An example of deformed mesh obtained after the solution 
of (12a) in the case of def A is shown in fig. |3(b)[ The deformed mesh is 
obtained by simply stretching each element according to the computed 
displacement and avoiding the insertion of any new vertex or further 
mesh refinement. 

Fig. shows the impact of channel size modification on the fluid 
velocity profiles (z-component) (no-ionic pressure is considered). Re¬ 
sults have been obtained by solving only ([^ in the deformed domain, 
with ID-cut taken along the geometrical symmetry axis. As expected, 
the average value of the velocity decreases with the increase of defor¬ 
mation, but in the region where the channel is restricted the electrolyte 
fluid velocity increases. 

Before describing the ion concentration profiles computed by the 
VE-PNP model in presence of mechanical forces it is important to 
look at the profiles of the driving forces acting on ion electrodiffusion. 
Figs. 12(a) and 12(b) illustrate the electrostatic potential and electric 
field profiles with and without deformation. Increasing the deformation 
the electric field strength in the channel center is quite reduced while 
the region close to SideA where the electric field is negative is extended 
towards the channel center. In all cases the electric field drop close 
to SideB corresponds to the almost flat profile in the electrostatic 
potential. 


The effects of ionic pressure on the fluid velocity are shown in fig. 13 


not only is the velocity not constant in the channel but also it undergoes 
a change of sign (inversion of electrolyte flow). From the cell-inside to 
the two third of the channel the velocity is directed towards the cell 
while in the final part of the channel it moves in the opposite direction. 











3D ION CHANNELS 


31 



Figure 11. Stokes equation in stationary conditions: 
ID-cnt along the geometrical symmetry axis of the z- 
component of the fluid velocity before and after mesh 
deformation. 




(a) (b) 

Figure 12. VE-PNP model in presence of mechanical 
deformation: ID-cut along the geometrical symmetry 
axis of the electrostatic potential (a) and of the electric 
held (b) at the final time T/ = 50ns. 


Moreover increasing deformation, a reduction of the velocity at the 
center of the channel is found due to the reduction of the electric field. 

Fig. M illustrates the K"*" and Na'*' concentration profiles. For both 
species the main variation with respect to the nnperturbed channel 
confignration is obtained in the center portion of the channel but with 
a remarkable difference. 
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Coordinate z [nm] 


Figure 13. VE-PNP model in presence of the mechani¬ 
cal deformation: ID-cut along the geometrical symmetry 
axis of the 2 :-component of the fluid velocity at the final 
time Tf = 50ns clearly highlight the effect of ionic pres¬ 
sure on the fluid. 


In the case of K'*', Pick’s diffusion is towards the end of the channel 
as the electric field and velocity in the case of normal size channel. 
This produces the K+ accumulation at the end of the channel: the 
Dirichlet boundary condition causes the reduction of the concentration 
at SideB. With the deformation, fluid velocity becomes negative and 
is directed towards the cell so that accumulation no longer occurs. 

In the case of Na"*", Pick’s diffusion is directed towards the cell while 
the velocity is towards the outside with an increased strength under 
deformation: this is the reason why the deformed channel has a less 
diffused profile. 



-'-1-' I ' I ' I ' I 

0 2 4 6 8 10 

Coordinate z [nm] 



- 1 —'- 1 —'- 1 —'- 1 —■— 
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(a) 


(b) 


Figure 14. VE-PNP model in presence of mechani¬ 
cal deformation: ID-cut along the geometrical symme¬ 
try axis of the ionic concentration at the final time 
Tf = 50ns: (a) K+, (b) Na+. 
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The 3D plots of fig. |15(a) show how the deformation is along the nor¬ 
mal direction of the channel axis and it increases with the deformation 
strength 7 in (44d) and (44e). As a consequence, channel diameter is 
further and further reduced. Fig. 15(b) shows how the potential profile 
of [T 2 (i 0 ] in a 3D view is uniform across the channel with no peaks 
induced by the possible accumulation of ionic charge at the channel 
surfaces. Fig. 16(a) shows the 3D view of the z-component of the 


fluid velocity: the profile across the cell is nonuniform also in the unde¬ 
formed case due to the adherence condition: the more the deformation, 
the more the nonuniformity of velocity distribution. Fig. |16(b) illus¬ 
trates the 3 D profile of Na+ ionic concentration. As in the case of 
electrolyte fluid, the increased channel deformation significantly affects 
the nonuniformity of the 3D distribution of the sodium ion inside the 
channel. This effect might be very important in regulating the elec¬ 
trostatic coupling between the mobile charge distribution in the chan¬ 
nel and the permanent charge that is usually trapped inside the lipid 
membrane bilayer surrounding the channel. This biophysical aspect of 
the problem has been investigated by Prof. W. Nonner and coworkers 
[7j and will be object of a future step of our research. 


m 



(a) (b) 

Figure 15. VE-PNP model in presence of mechanical 
deformation: 3D interior view of displacement (a) and 
electrostatic potential (b) at the final time T/ = 50ns. 

8. Conclusions and Future Perspectives 

In this article we have carried out a three-dimensional modeling and 
simulation, in three spatial dimensions and in time dependent condi¬ 
tions, of biological ion channels using a continuum-based approach. 
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(a) (b) 

Figure 16. VE-PNP model in presence of mechanical 
deformation: 3D interior view of the 2 :-component of fluid 
velocity (a) and of Na"*" concentration (b) at the flnal time 
Tf = 50ns. 


The proposed model is a multi-physics formulation that is capable to 
combine, to the best of our knowledge for the first time, ion electrod¬ 
iffusion, channel fluid motion, thermal self-heating and mechanical de¬ 
formation. The full self-consistency is intended for further study. Here, 
the thermal and mechanical effects are unidirectional. 

The resulting mathematical picture is a system of nonlinearly cou¬ 
pled partial differential equations in conservation form that describes 
the fundamental principles of ion charge and momentum balance, heat 
flow balance, electrolyte mass and momentum balance and mechanical 
balance in the nanochannel. 

A functional iteration that decouples the various differential sub¬ 
blocks of the system is introduced with the purpose of facilitating the 
computational implementation. Then, each resulting subproblem is 
discretized using the Galerkin Finite Element Method. 

The validation of the proposed computational model is carried out 
by simulating a realistic cylindrical voltage operated ion nanochannel 
where K+ and Na"*" ions are simultaneously flowing. Three sets of nu¬ 
merical experiments are performed and thoroughly discussed. We first 
investigate the coupling between electrochemical and fluid-dynamical 
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effects. Then, we enrich the modeling picture by investigating the influ¬ 
ence of a thermal gradient. Finally, we add a mechanical stress respon¬ 
sible of channel deformation and investigate its effect on the functional 
response of the channel. Results show that fluid and thermal fields have 
no influence in absence of mechanical deformation whereas ion distri¬ 
butions and channel functional response are significantly modified if 
mechanical stress is included in the model. These predictions agree 
with biophysical conjectures on the importance of protein conforma¬ 
tion in the modulation of channel electrochemical properties developed 
by Prof. W. Nonner and coworkers in [MllZl. 

Future extensions to overcome present model/methods limitations 
include, being not restricted to: 

• the coupling between stress and deformation fields and the PNP 
system to self-consistently calculate the deformation of channel 
wall; 

• a biophysically sound characterization of the heat production 
term Qheat the thermal block to account the effect of chemical 
reactions at the entrance of channel mouth between ions and, 
for example, externally supplied drugs; 

• the electrostatic coupling between the mobile charge distribu¬ 
tion in the channel and the permanent charge trapped inside 
the lipid membrane bilayer; 

• the extension and analysis of the EAFE formulation to account 
for electrolyte velocity in the thermo-electrical drift term in the 
discretization of the T-VE-PNP model. 
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